;Parallel Analysis in ViSta
;by Rubén Ledesma and Pedro Valero Mora
;e-mail rdledesma@gmail.com

;*PA-proto

(defproto pa-model-object-proto '(pca-or-fa method-value corcov-value cut-eigenvalues  obs-eigenvalues simulation percentil num-times nobs nvars) () analysis-plugin-object-proto)

(defmeth pa-model-object-proto :isnew 
  (menu-item-title tool-name model-prefix ok-data-types 
   data title name dialog ok-variable-types)
  (if (not (equal data *current-data*))(setcd data)
      
      (call-next-method 
             menu-item-title tool-name model-prefix ok-data-types 
             data title name dialog ok-variable-types)
      ))

;*SLOTS

(defmeth pa-model-object-proto :pca-or-fa
  (&optional (values nil set))
  (if set (setf (slot-value 'pca-or-fa) values))
  (slot-value 'pca-or-fa))
(defmeth pa-model-object-proto :obs-eigenvalues
  (&optional (values nil set))
  (if set (setf (slot-value 'obs-eigenvalues) values))
  (slot-value 'obs-eigenvalues))
(defmeth pa-model-object-proto :corcov-value
  (&optional (values nil set))
  (if set (setf (slot-value 'corcov-value) values)); 0=varianzas 1=correlaciones
  (slot-value 'corcov-value))
(defmeth pa-model-object-proto :method-value
  (&optional (values nil set))
  (if set (setf (slot-value 'method-value) values)); 0=simulate  1=permute
  (slot-value 'method-value))
(defmeth pa-model-object-proto :simulation
  (&optional (values nil set))
  (if set (setf (slot-value 'simulation) values))
  (slot-value 'simulation))
(defmeth pa-model-object-proto :cut-eigenvalues
  (&optional (values nil set))
  (if set (setf (slot-value 'cut-eigenvalues) values))
  (slot-value 'cut-eigenvalues))
(defmeth pa-model-object-proto :percentil
  (&optional (values nil set))
  (if set (setf (slot-value 'percentil) values))
  (slot-value 'percentil))
(defmeth pa-model-object-proto  :num-times
  (&optional (values nil set))
  (if set (setf (slot-value 'num-times) values))
  (slot-value 'num-times))
(defmeth pa-model-object-proto  :nobs
  (&optional (values nil set))
  (if set (setf (slot-value 'nobs) values))
  (slot-value 'nobs))
(defmeth pa-model-object-proto  :nvars
  (&optional (values nil set))
  (if set (setf (slot-value 'nvars) values))
  (slot-value 'nvars))

;*OPTIONS* 
(defmeth pa-model-object-proto :options ()

  (setf corcov (send choice-item-proto :new 
                            (list "Covariances"
                                  "Correlations")
                            :value 1))
  (setf method (send choice-item-proto :new 
                            (list "Simulate Normal Data"
                                  "Permute Data Matrix")
                            :value 0))
 (setf pca-or-fa (send choice-item-proto :new 
                            (list "Principal Component Analysis"
                                  "Principal Axis Factor Analysis")
                            :value 0))

  (setf text-percentil (send edit-text-item-proto :new ".95" :text-length 5))
  (setf text-times  (send edit-text-item-proto :new "100" :text-length 5))

  (setf OK (send modal-button-proto :new "  Ok"
                 :action
               #'(lambda ()
                   (let (
                         (dialog (send ok :dialog))
                         )
                     (send self :corcov-value (send corcov :value))
                     (send self :pca-or-fa (send pca-or-fa :value))
                     (send self :method-value (send method :value))
                     (send self :percentil (read-from-string (send text-percentil :text)))   
                     (send self :num-times (read-from-string (send text-times :text)))))))
                     
  (setf cancel (send modal-button-proto :new "  Cancel"
                     :action
               #'(lambda ()
                   (let (
                         (dialog (send cancel :dialog))
                         )
                     (send dialog :modal-dialog-return nil)))))

(setf vista-pa-dialog
      (send modal-dialog-proto :new 
            (list 
            ;(list "Compute from:" corcov)
             (list "  Model:    " pca-or-fa)
             (list "  Method:   " method)
             (list "  Cut-off percentil" text-percentil)
             (list "  Number of samples" text-times)
             (list ok cancel))
			:default-button ok :title "Options for Parallel Analysis"
            ))
(setf result (send vista-pa-dialog :modal-dialog))
             result)

;*ANALISIS

;;; Funciones de utilidad

 (defun my-eigen (m)
    (cond 
      ( 
        (= (send self :pca-or-fa) 1)
        (replace-diagonal m (- 1 (unicidad m)))           
        (setf my-eigen1 (eigenvalues  m))
        )
      (
        (= (send self :pca-or-fa) 0)
       (setf my-eigen (eigenvalues m))
       )
      ))

(defun unicidad (m)
  (/ (diagonal (inverse m))))

(defun replace-diagonal (m b)
  (let (
        (k (iseq (min (array-dimensions m))))
        )
    (mapcar #'(lambda (i)
                (setf (aref m i i) (elt b i))) k)
    )
  )   

(defun my-permute (data nvars nobs)
    (let* (
           (data data)
           (nvars nvars)
           (nobs nobs)           
(index (mapcar #'(lambda (x) (rank (UNIFORM-RAND x))) (repeat nobs nvars)))          
         )

      (setf  my-permute  (correlation-matrix  (transpose (make-array (list  nvars nobs ) :initial-contents
                                     (mapcar #'(lambda (column)
                                                 (select (select (column-list data) column) (order (select index column))))   (iseq nvars))))))
    )
  )

(defun n-permutations (n data nvars nobs)
    (mapcar #'(lambda (x) (my-permute data nvars nobs)) (iseq n))
    )


(defun sm0 (&optional (p 10) (m 100))
"Simula matrices de covarianzas utilizando la distribución de Wishart"
  (let* (
         (p p)
         (v (identity-matrix p))
         (m m)                
         (b-ji (make-array (list p p) :initial-element 0))
         (res-M nil)
         )
    (dotimes (i p)
             (dotimes (j p)
                      (when (> j i)
                            (setf (select b-ji i j) 
                                  (first (normal-rand 1))))
                      (when (= j i)
                            (setf (select b-ji j i) 
                                  (first (sqrt (chisq-rand 1 (+ 1 (- m (1+ j))))))))))
   (/ (cross-product b-ji) m)
    ))

(defun cov-to-corr (matriz)
     (let* (
            (matriz matriz)         
           (divisor (sqrt (outer-product (diagonal matriz) (diagonal matriz))))
            (res (/  matriz divisor)))
     res))

(defun sm (&optional (p 10) (m 100))
  (cov-to-corr (sm0 p m)))
  
(defmeth pa-model-object-proto :analysis () 
(let* (  
         (data (send current-data :active-data-matrix '(numeric)))
         (corcov-value (send self :corcov-value))
         (method-value (send self :method-value))
         (pca-or-fa (send self :pca-or-fa))
         (nvars (send current-data :active-nvar '(numeric)))
         (nobs (send current-data :active-nobs))
         (data (send current-data :active-data-matrix '(numeric)))
         (ntimes (send self :num-times))
         (percent (send self :percentil))
         (obs-eigenvalues (my-eigen (correlation-matrix data)))
         )

(defun n-my-eigen (ntimes nvars) 
      (make-array (list ntimes nvars) :initial-contents 
  (mapcar #'(lambda (x) (my-eigen x)) (n-permutations ntimes data nvars nobs))
            ))

(defun perc-my-eigen (p ntimes nvars)
(map-elements 'quantile 
                (column-list (n-my-eigen ntimes nvars)) p))

(defun compute-eigenvalues  (num-var nobs)
  (my-eigen (sm num-var nobs)))
       
 (defun MC-eigenvalues 
   (&key  (num-var 10) (corcov 0) (num-cases 100) (num-times 100) (num-eigenvalues 10))
   (apply 'bind-rows  
                     (mapcar #'(lambda (nt) 
                                 (select (compute-eigenvalues num-var num-cases) (iseq num-eigenvalues)))
                             (iseq num-times))))

(defun perc-MC-eigen 
  (&key  (num-var 10) (num-cases 100) (num-times 100) (num-eigenvalues 10) (percentil .95) (corcov  1))
  
  (map-elements 'quantile 
                (column-list (mc-eigenvalues 
                              :num-var num-var 
                              :num-cases num-cases 
                              :num-times num-times 
                              :num-eigenvalues num-eigenvalues
                              :corcov corcov))
                percentil))

(if (= method-value 0)
    (setf simulation 
(MC-eigenvalues :num-var nvars :num-cases nobs :num-times ntimes :num-eigenvalues nvars :corcov corcov-value)) (setf simulation (n-my-eigen ntimes nvars)))
    

(if (= method-value 0)
(setf cut-eigenvalues (perc-MC-eigen :num-var nvars :num-cases nobs :num-times ntimes :num-eigenvalues nvars :percentil percent :corcov corcov-value))

(setf cut-eigenvalues
      (perc-my-eigen percent ntimes nvars)))

  (send self :nobs nobs)
  (send self :nvars nvars)
  (send self :simulation simulation)
  (send self :cut-eigenvalues cut-eigenvalues)
  (send self :obs-eigenvalues obs-eigenvalues))
  )

;;;* REPORT 

(defmeth pa-model-object-proto :Report
  (&key (stream t) 
        (dialog nil))
;(if (not (eq current-object self)) (setcm self))
  (let* ((w nil)
         (var-labels (send self :variables))
         (nvars (send self :nvars))
         (nobs (send self :nobs))
         (Eig-Labels (mapcar #'(lambda (x) (format nil " Eigenvalue~a  " x))(+ 1 (iseq nvars))))
         (result (transpose (make-array (list 3 nvars) :initial-contents (list (send self :obs-eigenvalues) (map-elements 'mean 
(column-list (send self :simulation))) (send self :cut-eigenvalues) ))))
         )
          (setf w (report-header (send self :title "Parallel Analysis Report")))       
          (display-string
 (format nil "~%PARALLEL ANALYSIS REPORT~%by Ruben Ledesma & Pedro Valero Mora~2%MODEL: ~a~2%" (send self :name)) w)
          (display-string (format nil "VARIABLES: ~a~2%" var-labels ) w)


(display-string (format nil "Model: ~a~%" 
                                    (if (= 0 (send self  :pca-or-fa)) '("Principal Component of correlation matrix")
 '("Factor Analysis of correlation matrix"))) w)


(display-string (format nil "Method: ~a~%" 
                                    (if (= 0 (send self  :method-value)) '("Normal Data Simulation")
 '("Data Matrix Permutation"))) w)


    (display-string (format nil "Number of simulated samples: ~a~%" (send self :num-times) ) w)
    (display-string (format nil "Eigenvalues at percentile: ~a~%" (* 100 (send self :percentil))) w)

                (print-matrix-to-window result  w  :row-heading "" :row-labels Eig-Labels :column-heading "" :column-labels (list "Observed" "Mean" (format nil "Perc~a" (* 100 (send self :percentil)))) :decimals 5)
 (send w :fit-window-to-text)
 ))




;;*VISUALIZE


(defmeth pa-model-object-proto :visualize ()
  (let* (
         (data (send current-data :active-data-matrix '(numeric)))
         (nvars (send self :nvars))
         (nobs (send self :nvars))
         (ntimes (send self :num-times))
         (percent (send self :percentil))
         (cut-eigenvalues (send self :cut-eigenvalues))
         (simulation (send self :simulation))
         (obs-eigenvalues (send self :obs-eigenvalues))
         (Eig-Labels (mapcar #'(lambda (x) (format nil " Eigen.~a  " x))(+ 1 (iseq nvars))))

         (pa-plot (plot-points 
                 (iseq 1 nvars) cut-eigenvalues
                      :title "Parallel Plot"
                      :variable-labels '("" "Eigenvalue")
                      :show nil))

         (sim-box (boxplot (row-list (make-array (list nvars (+ 1 ntimes)) :initial-contents (combine (map-elements 'combine (send current-model :slot-value 'obs-eigenvalues) (column-list (send current-model :slot-value 'simulation)))))) :title "Simulation Box" :connect-points t :variable-labels nil :show nil))

         (sp (spreadplot (matrix '(1 2) 
                             (list  pa-plot sim-box))))
        )
    
    (send sim-box :SWITCH-CONNECT-POINTS)
    (send sim-box :new-menu "Sim-Box" 
              :items '(SHOWING-LABELS MOUSE DASH  
                            SYMBOL COLOR))
    (send sim-box :use-color t)
    (send sim-box :point-symbol (iseq ( + nvars (* ntimes nvars))) 'dot)
    (send sim-box :point-color (list 0) 'red)
    (send sim-box :point-state (iseq (+ (* ntimes nvars) nvars)) 'selected)
    (send sim-box :linked nil)
    (send sim-box :showing-labels nil)
    (send sim-box :add-lines (send sim-box :X) obs-eigenvalues :color 'red)
    (send sim-box :mouse-mode 'brushing)    
    (send sim-box :add-lines  (send sim-box :X) cut-eigenvalues :draw nil :color 'dark-green)


    (send pa-plot :add-lines    (list (iseq 1 nvars) cut-eigenvalues) :draw nil :color 'dark-green)
    (send pa-plot :point-color (iseq nvars) 'dark-green)
    (send pa-plot :add-lines    (list (iseq 1 nvars) (map-elements 'mean 
                                                                   (column-list simulation))) :draw nil :color 'grey)
    (send pa-plot :point-color (iseq nvars) 'dark-green)

    (mapcar #'(lambda (i)
                 (send pa-plot :point-label i
                       (cond
                         ((> i 0) 
                          (format nil "~5,4f" 
                                  (select cut-eigenvalues i) 
                                  ))
                         (t
                          (format nil "~5,4f" (select cut-eigenvalues i))))))
             (iseq (send pa-plot :num-points)))    
    (send pa-plot :range 1 nvars  1 :draw nil)
    (send pa-plot :plot-buttons :new-x nil :new-y nil :mouse-mode t)
    (send pa-plot :showing-labels t)
    (send pa-plot :mouse-mode 'brushing)
    (send pa-plot :menu nil)
    (send pa-plot :add-points (iseq 1 nvars) obs-eigenvalues)
    (send pa-plot :add-lines  (list (iseq 1 nvars) obs-eigenvalues) :draw nil :color 'red)
    (send pa-plot :point-color (iseq nvars (* 2 nvars)) 'red)

    (mapcar #'(lambda (i j)
              (send pa-plot :point-label i
                    (cond
                      ((> i 0) 
                       (format nil "~5,4f" 
                               (select obs-eigenvalues j) 
                               ))
                      (t
                       (format nil "~5,4f" (select obs-eigenvalues j))))))
            (iseq nvars (* 2 nvars)) (iseq nvars))
    (send pa-plot :adjust-to-data :draw nil)
    (send pa-plot :point-state (list 0 nvars) 'selected)
    (send sp :show-spreadplot))
  )
    
;;;;;;Create Data
(defmeth pa-model-object-proto :create-data 
  (&key (dialog nil)
        (simulation t)
        (mean-eigen nil)
        )

  (if (not (eq current-object self)) (setcm self))

  (let ((creator (send *desktop* :selected-icon))
        )
  (data (strcat "Eigen-" (send self :name))
   :created creator
   :creator-object self
   :title (strcat "Simulated eigenvalues for " (send self :title))
   :data  (combine (send self :simulation))
   :variables (mapcar #'(lambda (x) (format nil "Eign~a" x)) (iseq 1 (send self :nvars)))          
          :types (repeat "Numeric" (iseq (send self :nvars))))))

